temporary
2023-08-21
Clean data
Load data
dat <- readxl::read_xlsx("/Volumes/manhnd/Yi_roe_data.xlsx")
# filter variables for analyis
dat <- dat %>% dplyr::select("serialno", "country", "papm", "gender", "age", "education", "marital", "occupation", "live_area", "mobile", "internet")Transform data
# Remove the gender as others (18/2912)
dat1 <- dat %>%
filter(gender != 3 & gender != 4) %>%
mutate(marital = dplyr::recode(marital, `4` = 3, `5` = 3))
# Transform the valuables
dat1 <- dat1 %>%
mutate(
papm = factor(papm, ordered = TRUE),
gender = factor(gender, labels = c("Male", "Female")),
marital = factor(marital, labels = c("Never married", "Married/living with a partner", "Separated/widowed/divorced")),
education = factor(education, labels = c("No formal education", "Primary school", "Secondary school", "Diploma", "Degree")),
occupation = factor(occupation, labels = c("Employed", "Unemployed", "Student", "Retiree", "Homemaker")),
live_area = factor(live_area, labels = c("Urban", "Peri-urban", "Rural", "Slums")),
mobile = factor(mobile, labels = c("No", "Yes feature phone", "Yes smart phone")),
internet = factor(internet, labels = c("No", "Yes")),
country = factor(country)
)
summary(dat1)## serialno country papm gender age
## Min. : 1.0 Philippines:459 1:1596 Male :1365 Min. :18.00
## 1st Qu.: 726.2 Uganda :395 2: 291 Female:1529 1st Qu.:22.00
## Median :1453.5 Indonesia :368 3: 285 Median :30.00
## Mean :1455.6 Zimbabwe :292 4: 179 Mean :34.95
## 3rd Qu.:2185.8 Cameroon :289 5: 543 3rd Qu.:44.00
## Max. :2912.0 Bangladesh :280 Max. :90.00
## (Other) :811
## education marital
## No formal education: 173 Never married :1221
## Primary school : 476 Married/living with a partner:1395
## Secondary school :1182 Separated/widowed/divorced : 278
## Diploma : 405
## Degree : 658
##
##
## occupation live_area mobile internet
## Employed :1294 Urban :1066 No : 257 No : 839
## Unemployed: 566 Peri-urban: 374 Yes feature phone: 626 Yes:2055
## Student : 681 Rural :1281 Yes smart phone :2011
## Retiree : 65 Slums : 173
## Homemaker : 288
##
## str(dat1)## tibble [2,894 × 11] (S3: tbl_df/tbl/data.frame)
## $ serialno : num [1:2894] 1 2 3 4 5 6 7 8 9 10 ...
## $ country : Factor w/ 9 levels "Bangladesh","Cameroon",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ papm : Ord.factor w/ 5 levels "1"<"2"<"3"<"4"<..: 1 1 1 1 1 1 1 1 1 1 ...
## $ gender : Factor w/ 2 levels "Male","Female": 2 2 2 2 2 1 1 2 1 2 ...
## $ age : num [1:2894] 29 36 52 27 31 31 70 50 49 45 ...
## $ education : Factor w/ 5 levels "No formal education",..: 2 1 1 2 2 2 1 1 1 1 ...
## $ marital : Factor w/ 3 levels "Never married",..: 2 2 3 2 2 2 2 2 2 2 ...
## $ occupation: Factor w/ 5 levels "Employed","Unemployed",..: 5 5 1 5 5 1 2 5 1 5 ...
## $ live_area : Factor w/ 4 levels "Urban","Peri-urban",..: 3 3 3 3 3 3 3 3 3 3 ...
## $ mobile : Factor w/ 3 levels "No","Yes feature phone",..: 2 2 2 2 2 2 2 1 2 2 ...
## $ internet : Factor w/ 2 levels "No","Yes": 1 1 1 1 1 1 1 1 1 1 ...Visualize some variables
papm versus countries
dat2 <- dat1 %>%
select("country", "papm") %>%
group_by(country, papm) %>%
summarise(n = n()) %>%
mutate(proportion_papm = n / sum(n)) %>%
ungroup()## `summarise()` has grouped output by 'country'. You can override using the
## `.groups` argument.ggplot(data = dat2, aes(x = country, y = proportion_papm)) +
geom_col(width = 0.5) +
facet_grid(~papm) +
coord_flip()
Education versus countries
dat2 <- dat %>%
select("country", "education") %>%
group_by(country, education) %>%
summarise(n = n()) %>%
mutate(proportion_education = n / sum(n)) %>%
ungroup()## `summarise()` has grouped output by 'country'. You can override using the
## `.groups` argument.ggplot(data = dat2, aes(x = country, y = proportion_education)) +
geom_col(width = 0.5) +
facet_grid(~education) +
coord_flip()
MASS package: polr modell
Fit the modell without country interaction
m1 <- polr(formula = papm ~ country + gender + age + education + marital + occupation + live_area + mobile + internet, data = dat1, Hess = TRUE)
summary(m1)## Call:
## polr(formula = papm ~ country + gender + age + education + marital +
## occupation + live_area + mobile + internet, data = dat1,
## Hess = TRUE)
##
## Coefficients:
## Value Std. Error t value
## countryCameroon 1.030800 0.304414 3.38618
## countryIndia 3.224494 0.272151 11.84819
## countryIndonesia 1.723045 0.279535 6.16397
## countryKenya 2.155655 0.310191 6.94944
## countryNepal 3.933541 0.282818 13.90837
## countryPhilippines 1.061534 0.286374 3.70681
## countryUganda 4.131868 0.272311 15.17337
## countryZimbabwe 4.363371 0.287158 15.19501
## genderFemale -0.005772 0.085817 -0.06726
## age -0.001336 0.003866 -0.34557
## educationPrimary school 0.480407 0.190301 2.52446
## educationSecondary school 0.858222 0.195790 4.38338
## educationDiploma 1.083491 0.222924 4.86036
## educationDegree 1.204571 0.220266 5.46872
## maritalMarried/living with a partner 0.243062 0.122772 1.97978
## maritalSeparated/widowed/divorced 0.602561 0.181871 3.31312
## occupationUnemployed -0.183360 0.118656 -1.54530
## occupationStudent -0.156333 0.134595 -1.16151
## occupationRetiree -0.556455 0.284433 -1.95637
## occupationHomemaker -0.404096 0.151572 -2.66603
## live_areaPeri-urban -0.019113 0.135460 -0.14110
## live_areaRural 0.056872 0.116823 0.48682
## live_areaSlums 0.341699 0.225495 1.51533
## mobileYes feature phone 0.576591 0.161966 3.55994
## mobileYes smart phone 0.217389 0.212906 1.02106
## internetYes 0.491155 0.168225 2.91964
##
## Intercepts:
## Value Std. Error t value
## 1|2 4.2778 0.3845 11.1264
## 2|3 4.9192 0.3873 12.7018
## 3|4 5.5729 0.3902 14.2807
## 4|5 6.0354 0.3921 15.3938
##
## Residual Deviance: 6159.317
## AIC: 6219.317Fit the model with country interaction
# First fit the model with interaction between country and all of the independent variables but the model cannot run. Then I removed the interaction between live_area and country, model can run as below
m2 <- polr(formula = papm ~ country + gender + age + education + marital + occupation + live_area + mobile + internet + gender:country + age:country + marital:country + occupation:country + mobile:country + internet:country + education:country, data = dat1, Hess = TRUE)## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurredsummary(m2)## Call:
## polr(formula = papm ~ country + gender + age + education + marital +
## occupation + live_area + mobile + internet + gender:country +
## age:country + marital:country + occupation:country + mobile:country +
## internet:country + education:country, data = dat1, Hess = TRUE)
##
## Coefficients:
## Value Std. Error
## countryCameroon 14.609312 1.400e+00
## countryIndia 24.913158 1.278e+00
## countryIndonesia -3.055593 1.099e+00
## countryKenya 26.698692 1.697e+00
## countryNepal 27.589432 1.359e+00
## countryPhilippines 23.911632 1.497e+00
## countryUganda 26.874774 1.138e+00
## countryZimbabwe 16.276852 1.281e+00
## genderFemale 0.354487 5.578e-01
## age 0.077100 3.809e-02
## educationPrimary school -4.080201 3.335e-01
## educationSecondary school 9.866479 7.889e-01
## educationDiploma 11.217361 6.772e-01
## educationDegree 11.537265 5.694e-01
## maritalMarried/living with a partner -0.857760 1.184e+00
## maritalSeparated/widowed/divorced -2.453037 2.089e-01
## occupationUnemployed -16.542871 1.461e-01
## occupationStudent 0.048204 1.308e+00
## occupationRetiree -18.167315 3.829e-01
## occupationHomemaker -12.483776 2.274e-01
## live_areaPeri-urban -0.090668 1.443e-01
## live_areaRural -0.001955 1.294e-01
## live_areaSlums 0.249086 2.547e-01
## mobileYes feature phone 11.977607 8.830e-01
## mobileYes smart phone 11.867155 1.091e+00
## internetYes 0.097749 1.624e+00
## countryCameroon:genderFemale -0.475020 6.699e-01
## countryIndia:genderFemale 0.074489 6.111e-01
## countryIndonesia:genderFemale -0.611408 6.207e-01
## countryKenya:genderFemale -0.179755 6.525e-01
## countryNepal:genderFemale -0.674883 6.148e-01
## countryPhilippines:genderFemale -0.091360 6.396e-01
## countryUganda:genderFemale -0.477477 5.858e-01
## countryZimbabwe:genderFemale -0.918931 6.375e-01
## countryCameroon:age -0.074785 4.493e-02
## countryIndia:age -0.063382 3.976e-02
## countryIndonesia:age -0.033517 4.247e-02
## countryKenya:age -0.031155 4.203e-02
## countryNepal:age -0.064597 4.016e-02
## countryPhilippines:age -0.080832 4.119e-02
## countryUganda:age -0.073202 3.865e-02
## countryZimbabwe:age -0.122026 4.188e-02
## countryCameroon:maritalMarried/living with a partner 1.420520 1.273e+00
## countryIndia:maritalMarried/living with a partner 0.976123 1.236e+00
## countryIndonesia:maritalMarried/living with a partner 0.718878 1.232e+00
## countryKenya:maritalMarried/living with a partner 0.783927 1.265e+00
## countryNepal:maritalMarried/living with a partner 0.381139 1.243e+00
## countryPhilippines:maritalMarried/living with a partner 1.822441 1.282e+00
## countryUganda:maritalMarried/living with a partner 0.919528 1.221e+00
## countryZimbabwe:maritalMarried/living with a partner 2.381148 1.259e+00
## countryCameroon:maritalSeparated/widowed/divorced -12.787097 8.655e-08
## countryIndia:maritalSeparated/widowed/divorced 2.628711 6.726e-01
## countryIndonesia:maritalSeparated/widowed/divorced 0.975907 6.736e-01
## countryKenya:maritalSeparated/widowed/divorced 2.709081 5.056e-01
## countryNepal:maritalSeparated/widowed/divorced 3.123797 6.262e-01
## countryPhilippines:maritalSeparated/widowed/divorced 3.826500 6.370e-01
## countryUganda:maritalSeparated/widowed/divorced 2.576574 4.045e-01
## countryZimbabwe:maritalSeparated/widowed/divorced 4.439632 5.209e-01
## countryCameroon:occupationUnemployed 16.766464 4.480e-01
## countryIndia:occupationUnemployed 16.584719 4.851e-01
## countryIndonesia:occupationUnemployed 15.559812 4.244e-01
## countryKenya:occupationUnemployed 16.774798 4.632e-01
## countryNepal:occupationUnemployed 14.736578 5.290e-01
## countryPhilippines:occupationUnemployed 16.481324 4.856e-01
## countryUganda:occupationUnemployed 16.816579 2.401e-01
## countryZimbabwe:occupationUnemployed 15.691493 3.303e-01
## countryCameroon:occupationStudent -0.162487 1.393e+00
## countryIndia:occupationStudent 0.259397 1.358e+00
## countryIndonesia:occupationStudent -0.397042 1.361e+00
## countryKenya:occupationStudent -0.375236 1.443e+00
## countryNepal:occupationStudent -0.461873 1.359e+00
## countryPhilippines:occupationStudent 0.709372 1.412e+00
## countryUganda:occupationStudent 0.011918 1.360e+00
## countryZimbabwe:occupationStudent 0.052282 1.394e+00
## countryCameroon:occupationRetiree 0.304758 1.344e-08
## countryIndia:occupationRetiree 17.765439 6.078e-01
## countryIndonesia:occupationRetiree 1.689434 1.391e-07
## countryKenya:occupationRetiree 18.059350 1.465e+00
## countryNepal:occupationRetiree 16.820388 5.936e-01
## countryPhilippines:occupationRetiree 18.007567 8.731e-01
## countryUganda:occupationRetiree 17.472933 7.974e-01
## countryZimbabwe:occupationRetiree 36.326893 6.540e-08
## countryCameroon:occupationHomemaker 11.750979 9.985e-01
## countryIndia:occupationHomemaker 11.743691 3.820e-01
## countryIndonesia:occupationHomemaker 11.075414 6.715e-01
## countryKenya:occupationHomemaker 13.993172 7.158e-01
## countryNepal:occupationHomemaker 11.784830 4.632e-01
## countryPhilippines:occupationHomemaker 12.831394 6.032e-01
## countryUganda:occupationHomemaker 12.899331 3.283e-01
## countryZimbabwe:occupationHomemaker 10.611018 9.371e-01
## countryCameroon:mobileYes feature phone -12.848847 1.267e+00
## countryIndia:mobileYes feature phone -11.411442 9.488e-01
## countryIndonesia:mobileYes feature phone -14.822680 2.206e+00
## countryKenya:mobileYes feature phone -11.090439 1.084e+00
## countryNepal:mobileYes feature phone -12.424780 1.400e+00
## countryPhilippines:mobileYes feature phone -11.984637 1.099e+00
## countryUganda:mobileYes feature phone -11.399280 9.077e-01
## countryZimbabwe:mobileYes feature phone -8.440743 1.168e+00
## countryCameroon:mobileYes smart phone -11.923094 1.571e+00
## countryIndia:mobileYes smart phone -11.281882 1.357e+00
## countryIndonesia:mobileYes smart phone -13.735593 2.200e+00
## countryKenya:mobileYes smart phone -10.639113 1.271e+00
## countryNepal:mobileYes smart phone -12.697334 1.963e+00
## countryPhilippines:mobileYes smart phone -13.246191 1.279e+00
## countryUganda:mobileYes smart phone -11.300390 1.198e+00
## countryZimbabwe:mobileYes smart phone -10.451132 1.268e+00
## countryCameroon:internetYes -0.142093 1.898e+00
## countryIndia:internetYes 0.343617 1.800e+00
## countryIndonesia:internetYes 19.075766 1.099e+00
## countryKenya:internetYes -0.298728 1.683e+00
## countryNepal:internetYes 1.112136 2.262e+00
## countryPhilippines:internetYes 2.030549 1.746e+00
## countryUganda:internetYes 0.590445 1.694e+00
## countryZimbabwe:internetYes -0.705238 1.672e+00
## countryCameroon:educationPrimary school 13.077706 6.222e-01
## countryIndia:educationPrimary school 4.762476 5.766e-01
## countryIndonesia:educationPrimary school 16.204295 7.376e-01
## countryKenya:educationPrimary school 0.185179 1.174e+00
## countryNepal:educationPrimary school 3.920423 5.980e-01
## countryPhilippines:educationPrimary school 3.667566 8.605e-01
## countryUganda:educationPrimary school 4.581924 3.971e-01
## countryZimbabwe:educationPrimary school 17.735081 7.247e-01
## countryCameroon:educationSecondary school 0.958803 8.201e-01
## countryIndia:educationSecondary school -8.921375 9.204e-01
## countryIndonesia:educationSecondary school 1.257021 7.397e-01
## countryKenya:educationSecondary school -12.508199 1.331e+00
## countryNepal:educationSecondary school -9.343525 9.348e-01
## countryPhilippines:educationSecondary school -9.735912 1.081e+00
## countryUganda:educationSecondary school -9.409909 8.247e-01
## countryZimbabwe:educationSecondary school 3.027823 7.846e-01
## countryCameroon:educationDiploma 0.244220 7.472e-01
## countryIndia:educationDiploma -9.866073 8.812e-01
## countryIndonesia:educationDiploma 0.688810 6.993e-01
## countryKenya:educationDiploma -14.272035 1.290e+00
## countryNepal:educationDiploma -9.862474 8.626e-01
## countryPhilippines:educationDiploma -11.230641 1.041e+00
## countryUganda:educationDiploma -10.547329 7.987e-01
## countryZimbabwe:educationDiploma 1.256645 7.247e-01
## countryCameroon:educationDegree -0.181528 6.094e-01
## countryIndia:educationDegree -9.934391 7.724e-01
## countryIndonesia:educationDegree 0.151529 5.295e-01
## countryKenya:educationDegree -13.209072 1.288e+00
## countryNepal:educationDegree -10.477752 7.820e-01
## countryPhilippines:educationDegree -11.690755 9.826e-01
## countryUganda:educationDegree -11.004074 1.124e+00
## countryZimbabwe:educationDegree 0.496835 5.932e-01
## t value
## countryCameroon 1.043e+01
## countryIndia 1.950e+01
## countryIndonesia -2.779e+00
## countryKenya 1.574e+01
## countryNepal 2.030e+01
## countryPhilippines 1.597e+01
## countryUganda 2.362e+01
## countryZimbabwe 1.270e+01
## genderFemale 6.355e-01
## age 2.024e+00
## educationPrimary school -1.224e+01
## educationSecondary school 1.251e+01
## educationDiploma 1.656e+01
## educationDegree 2.026e+01
## maritalMarried/living with a partner -7.245e-01
## maritalSeparated/widowed/divorced -1.174e+01
## occupationUnemployed -1.132e+02
## occupationStudent 3.684e-02
## occupationRetiree -4.744e+01
## occupationHomemaker -5.490e+01
## live_areaPeri-urban -6.285e-01
## live_areaRural -1.511e-02
## live_areaSlums 9.779e-01
## mobileYes feature phone 1.356e+01
## mobileYes smart phone 1.088e+01
## internetYes 6.018e-02
## countryCameroon:genderFemale -7.091e-01
## countryIndia:genderFemale 1.219e-01
## countryIndonesia:genderFemale -9.851e-01
## countryKenya:genderFemale -2.755e-01
## countryNepal:genderFemale -1.098e+00
## countryPhilippines:genderFemale -1.428e-01
## countryUganda:genderFemale -8.151e-01
## countryZimbabwe:genderFemale -1.442e+00
## countryCameroon:age -1.664e+00
## countryIndia:age -1.594e+00
## countryIndonesia:age -7.893e-01
## countryKenya:age -7.413e-01
## countryNepal:age -1.609e+00
## countryPhilippines:age -1.962e+00
## countryUganda:age -1.894e+00
## countryZimbabwe:age -2.914e+00
## countryCameroon:maritalMarried/living with a partner 1.116e+00
## countryIndia:maritalMarried/living with a partner 7.898e-01
## countryIndonesia:maritalMarried/living with a partner 5.834e-01
## countryKenya:maritalMarried/living with a partner 6.198e-01
## countryNepal:maritalMarried/living with a partner 3.066e-01
## countryPhilippines:maritalMarried/living with a partner 1.421e+00
## countryUganda:maritalMarried/living with a partner 7.534e-01
## countryZimbabwe:maritalMarried/living with a partner 1.892e+00
## countryCameroon:maritalSeparated/widowed/divorced -1.477e+08
## countryIndia:maritalSeparated/widowed/divorced 3.908e+00
## countryIndonesia:maritalSeparated/widowed/divorced 1.449e+00
## countryKenya:maritalSeparated/widowed/divorced 5.359e+00
## countryNepal:maritalSeparated/widowed/divorced 4.989e+00
## countryPhilippines:maritalSeparated/widowed/divorced 6.007e+00
## countryUganda:maritalSeparated/widowed/divorced 6.369e+00
## countryZimbabwe:maritalSeparated/widowed/divorced 8.523e+00
## countryCameroon:occupationUnemployed 3.743e+01
## countryIndia:occupationUnemployed 3.419e+01
## countryIndonesia:occupationUnemployed 3.666e+01
## countryKenya:occupationUnemployed 3.622e+01
## countryNepal:occupationUnemployed 2.786e+01
## countryPhilippines:occupationUnemployed 3.394e+01
## countryUganda:occupationUnemployed 7.003e+01
## countryZimbabwe:occupationUnemployed 4.750e+01
## countryCameroon:occupationStudent -1.166e-01
## countryIndia:occupationStudent 1.910e-01
## countryIndonesia:occupationStudent -2.918e-01
## countryKenya:occupationStudent -2.601e-01
## countryNepal:occupationStudent -3.400e-01
## countryPhilippines:occupationStudent 5.024e-01
## countryUganda:occupationStudent 8.762e-03
## countryZimbabwe:occupationStudent 3.751e-02
## countryCameroon:occupationRetiree 2.267e+07
## countryIndia:occupationRetiree 2.923e+01
## countryIndonesia:occupationRetiree 1.215e+07
## countryKenya:occupationRetiree 1.233e+01
## countryNepal:occupationRetiree 2.834e+01
## countryPhilippines:occupationRetiree 2.063e+01
## countryUganda:occupationRetiree 2.191e+01
## countryZimbabwe:occupationRetiree 5.555e+08
## countryCameroon:occupationHomemaker 1.177e+01
## countryIndia:occupationHomemaker 3.075e+01
## countryIndonesia:occupationHomemaker 1.649e+01
## countryKenya:occupationHomemaker 1.955e+01
## countryNepal:occupationHomemaker 2.544e+01
## countryPhilippines:occupationHomemaker 2.127e+01
## countryUganda:occupationHomemaker 3.929e+01
## countryZimbabwe:occupationHomemaker 1.132e+01
## countryCameroon:mobileYes feature phone -1.014e+01
## countryIndia:mobileYes feature phone -1.203e+01
## countryIndonesia:mobileYes feature phone -6.720e+00
## countryKenya:mobileYes feature phone -1.023e+01
## countryNepal:mobileYes feature phone -8.874e+00
## countryPhilippines:mobileYes feature phone -1.090e+01
## countryUganda:mobileYes feature phone -1.256e+01
## countryZimbabwe:mobileYes feature phone -7.224e+00
## countryCameroon:mobileYes smart phone -7.591e+00
## countryIndia:mobileYes smart phone -8.314e+00
## countryIndonesia:mobileYes smart phone -6.244e+00
## countryKenya:mobileYes smart phone -8.372e+00
## countryNepal:mobileYes smart phone -6.467e+00
## countryPhilippines:mobileYes smart phone -1.035e+01
## countryUganda:mobileYes smart phone -9.429e+00
## countryZimbabwe:mobileYes smart phone -8.240e+00
## countryCameroon:internetYes -7.487e-02
## countryIndia:internetYes 1.909e-01
## countryIndonesia:internetYes 1.735e+01
## countryKenya:internetYes -1.775e-01
## countryNepal:internetYes 4.916e-01
## countryPhilippines:internetYes 1.163e+00
## countryUganda:internetYes 3.485e-01
## countryZimbabwe:internetYes -4.219e-01
## countryCameroon:educationPrimary school 2.102e+01
## countryIndia:educationPrimary school 8.260e+00
## countryIndonesia:educationPrimary school 2.197e+01
## countryKenya:educationPrimary school 1.578e-01
## countryNepal:educationPrimary school 6.556e+00
## countryPhilippines:educationPrimary school 4.262e+00
## countryUganda:educationPrimary school 1.154e+01
## countryZimbabwe:educationPrimary school 2.447e+01
## countryCameroon:educationSecondary school 1.169e+00
## countryIndia:educationSecondary school -9.693e+00
## countryIndonesia:educationSecondary school 1.699e+00
## countryKenya:educationSecondary school -9.395e+00
## countryNepal:educationSecondary school -9.995e+00
## countryPhilippines:educationSecondary school -9.008e+00
## countryUganda:educationSecondary school -1.141e+01
## countryZimbabwe:educationSecondary school 3.859e+00
## countryCameroon:educationDiploma 3.269e-01
## countryIndia:educationDiploma -1.120e+01
## countryIndonesia:educationDiploma 9.850e-01
## countryKenya:educationDiploma -1.107e+01
## countryNepal:educationDiploma -1.143e+01
## countryPhilippines:educationDiploma -1.079e+01
## countryUganda:educationDiploma -1.321e+01
## countryZimbabwe:educationDiploma 1.734e+00
## countryCameroon:educationDegree -2.979e-01
## countryIndia:educationDegree -1.286e+01
## countryIndonesia:educationDegree 2.862e-01
## countryKenya:educationDegree -1.025e+01
## countryNepal:educationDegree -1.340e+01
## countryPhilippines:educationDegree -1.190e+01
## countryUganda:educationDegree -9.794e+00
## countryZimbabwe:educationDegree 8.376e-01
##
## Intercepts:
## Value Std. Error t value
## 1|2 2.706380e+01 1.076900e+00 2.513160e+01
## 2|3 2.775100e+01 1.077600e+00 2.575330e+01
## 3|4 2.845470e+01 1.078200e+00 2.638990e+01
## 4|5 2.895350e+01 1.078700e+00 2.684100e+01
##
## Residual Deviance: 5866.142
## AIC: 6166.142# compare with model 1
anova(m1, m2, test = "Chisq")## Likelihood ratio tests of ordinal regression models
##
## Response: papm
## Model
## 1 country + gender + age + education + marital + occupation + live_area + mobile + internet
## 2 country + gender + age + education + marital + occupation + live_area + mobile + internet + gender:country + age:country + marital:country + occupation:country + mobile:country + internet:country + education:country
## Resid. df Resid. Dev Test Df LR stat. Pr(Chi)
## 1 2864 6159.317
## 2 2744 5866.142 1 vs 2 120 293.1756 1.110223e-16Anova(m2)## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred## Analysis of Deviance Table (Type II tests)
##
## Response: papm
## LR Chisq Df Pr(>Chisq)
## country 1052.99 8 < 2.2e-16 ***
## gender 0.41 1 0.5230524
## age 3.39 1 0.0654129 .
## education 28.66 4 9.161e-06 ***
## marital 3.83 2 0.1470023
## occupation 7.84 4 0.0975951 .
## live_area 1.64 3 0.6506180
## mobile 15.12 2 0.0005211 ***
## internet 2.95 1 0.0860821 .
## country:gender 10.36 8 0.2405672
## country:age 22.71 8 0.0037536 **
## country:marital 31.38 16 0.0120169 *
## country:occupation 59.74 32 0.0020828 **
## country:mobile 41.60 16 0.0004521 ***
## country:internet 28.18 8 0.0004409 ***
## country:education 70.05 32 0.0001168 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# m2 is beter than m1.Testing the proportional odds assumption
https://peopleanalytics-regression-book.org/ord-reg.html:
dat_test1 <- dat1 %>% mutate(
test1 = ifelse(papm == "1", 0, 1),
test2 = ifelse(papm == "1" | papm == "2", 0, 1),
test3 = ifelse(papm == "1" | papm == "2" | papm == "3", 0, 1),
test4 = ifelse(papm == "1" | papm == "2" | papm == "3" | papm == "4", 0, 1)
)
table(dat_test1$papm)##
## 1 2 3 4 5
## 1596 291 285 179 543m1 <- glm(test1 ~ gender + age + education + marital + occupation + live_area + mobile + internet,
data = dat_test1 %>% filter(country == "Bangladesh"),
family = "binomial"
)## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurredm2 <- glm(test2 ~ gender + age + education + marital + occupation + live_area + mobile + internet,
data = dat_test1 %>% filter(country == "Bangladesh"),
family = "binomial"
)## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurredm3 <- glm(test3 ~ gender + age + education + marital + occupation + live_area + mobile + internet,
data = dat_test1 %>% filter(country == "Bangladesh"),
family = "binomial"
)## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurredm4 <- glm(test4 ~ gender + age + education + marital + occupation + live_area + mobile + internet,
data = dat_test1 %>% filter(country == "Bangladesh"),
family = "binomial"
)## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurredcoefficient_comparison <- data.frame(
test1 = summary(m1)$coefficients[, "Estimate"],
test2 = summary(m2)$coefficients[, "Estimate"],
test3 = summary(m3)$coefficients[, "Estimate"],
test4 = summary(m4)$coefficients[, "Estimate"],
diff1_2 = summary(m2)$coefficients[, "Estimate"] - summary(m1)$coefficients[, "Estimate"],
diff1_3 = summary(m3)$coefficients[, "Estimate"] - summary(m1)$coefficients[, "Estimate"],
diff1_4 = summary(m4)$coefficients[, "Estimate"] - summary(m1)$coefficients[, "Estimate"]
)
coefficient_comparison## test1 test2 test3
## (Intercept) -38.43425101 -38.43425101 -40.93666714
## genderFemale 0.38700242 0.38700242 0.09338424
## age 0.08512152 0.08512152 0.08291879
## educationPrimary school 1.53180401 1.53180401 1.49493880
## educationSecondary school 17.77258036 17.77258036 18.95926562
## educationDiploma 19.23192113 19.23192113 19.17753807
## educationDegree 19.54783604 19.54783604 20.62539567
## maritalMarried/living with a partner -0.65031291 -0.65031291 -0.56711481
## maritalSeparated/widowed/divorced -2.11973287 -2.11973287 -1.11113159
## occupationUnemployed -18.58084035 -18.58084035 -19.58069143
## occupationStudent 0.11277406 0.11277406 -0.67586541
## occupationRetiree -20.44333003 -20.44333003 -20.16806042
## occupationHomemaker -16.21747376 -16.21747376 -16.88829039
## live_areaPeri-urban -0.45811434 -0.45811434 1.00093169
## live_areaRural 0.77632174 0.77632174 2.28880765
## mobileYes feature phone 14.16567476 14.16567476 14.22825219
## mobileYes smart phone 14.11045566 14.11045566 14.20706684
## internetYes 0.28136377 0.28136377 0.38110277
## test4 diff1_2 diff1_3
## (Intercept) -61.18243089 0 -2.502416133
## genderFemale -0.11672626 0 -0.293618176
## age 0.03502682 0 -0.002202736
## educationPrimary school 1.07945192 0 -0.036865204
## educationSecondary school 19.16571382 0 1.186685264
## educationDiploma -2.44021033 0 -0.054383058
## educationDegree 20.91176452 0 1.077559635
## maritalMarried/living with a partner 17.96932349 0 0.083198108
## maritalSeparated/widowed/divorced 18.74805329 0 1.008601275
## occupationUnemployed -16.15798504 0 -0.999851084
## occupationStudent -9.07947685 0 -0.788639469
## occupationRetiree -1.73019043 0 0.275269612
## occupationHomemaker -17.96936153 0 -0.670816628
## live_areaPeri-urban 18.80811464 0 1.459046033
## live_areaRural 20.42660607 0 1.512485918
## mobileYes feature phone -0.59054383 0 0.062577434
## mobileYes smart phone 5.33094271 0 0.096611175
## internetYes -5.93206007 0 0.099739006
## diff1_4
## (Intercept) -22.74817988
## genderFemale -0.50372868
## age -0.05009471
## educationPrimary school -0.45235209
## educationSecondary school 1.39313347
## educationDiploma -21.67213146
## educationDegree 1.36392848
## maritalMarried/living with a partner 18.61963640
## maritalSeparated/widowed/divorced 20.86778615
## occupationUnemployed 2.42285531
## occupationStudent -9.19225091
## occupationRetiree 18.71313961
## occupationHomemaker -1.75188777
## live_areaPeri-urban 19.26622898
## live_areaRural 19.65028434
## mobileYes feature phone -14.75621859
## mobileYes smart phone -8.77951296
## internetYes -6.21342384We could see there is a large difference coefficients. The assumption of proportional odds is violated. We need to use different methods. The book recommend some other methods explain more details in Agresti (2010))
- Baseline logistic model. This model is the same as the multinomial regression model covered in the previous chapter, using the lowest ordinal value as the reference.
- Adjacent-category logistic model. This model compares each level of the ordinal variable to the next highest level, and it is a constrained version of the baseline logistic model. The brglm2 package in R offers a function bracl() for calculating an adjacent category logistic model.
- Continuation-ratio logistic model. This model compares each level of the ordinal variable to all lower levels. This can be modeled using binary logistic regression techniques, but new variables need to be constructed from the data set to allow this. The R package rms has a function cr.setup() which is a utility for preparing an outcome variable for a continuation ratio model.
I still fit the model to individual country to see where the issues occur. When the absolute Value of the predictor is greater than 10, I make table for that independent variable and the outcome for that country. ## Calculate the confident interval I use Wald test to obtain the CI because the likelihood ratio test doesn’t work
countryList <- c("Bangladesh", "Cameroon", "India", "Indonesia", "Kenya", "Nepal", "Philippines", "Uganda", "Zimbabwe")
fun <- function(data, countryList) {
allData <- list()
for (ctry in countryList) {
m <- polr(formula = papm ~ gender + age + education + marital + occupation + live_area + mobile + internet, data = dat1 %>% filter(country == ctry), Hess = TRUE)
ci <- confint.default(m)
tmp <- as.data.frame(round(exp(cbind(OR = coef(m), ci)), 2)) %>%
mutate(
country = ctry,
predictor = row.names(m)
)
allData[[ctry]] <- tmp
}
combinedData <- do.call(rbind, allData)
return(combinedData)
}
tmp <- fun(dat1, countryList)## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred## Warning in polr(formula = papm ~ gender + age + education + marital + occupation
## + : design appears to be rank-deficient, so dropping some coefs## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred## Warning in polr(formula = papm ~ gender + age + education + marital + occupation
## + : design appears to be rank-deficient, so dropping some coefs
## Warning in polr(formula = papm ~ gender + age + education + marital + occupation
## + : design appears to be rank-deficient, so dropping some coefsVisualize the result
Clean the result
tmp1 <- tmp %>%
mutate(independant_var = row.names(tmp)) %>%
separate(independant_var, c("count", "var"), sep = "\\.") %>%
mutate(var = str_replace_all(var, c(
"gender" = "gender.",
"age" = "age.Age",
"education" = "education.",
"marital" = "marital.",
"occupation" = "occupation.",
"live_area" = "live_area.",
"mobile" = "mobile.",
"internet" = "internet."
))) %>%
separate(var, c("var", "level"), "\\.") %>%
mutate(
OR = ifelse(OR >= 1000, 0, OR),
`2.5 %` = ifelse(`2.5 %` >= 1000, 0, `2.5 %`),
`97.5 %` = ifelse(`97.5 %` >= 1000, 0, `97.5 %`)
) %>%
mutate(level = as.factor(level))Gender
tmp2 <- tmp1 %>%
filter(var == "gender")
p <- ggplot(data = tmp2, aes(x = level, y = OR, colour = country)) +
geom_line(group = 1, position = position_dodge(0.3)) +
geom_errorbar(aes(ymin = `2.5 %`, ymax = `97.5 %`), position = position_dodge(0.3), width = 2) +
geom_point(size = 3, position = position_dodge(0.3)) +
scale_y_continuous(trans = "log10") +
scale_fill_brewer(palette = "Set1")+
labs(title = "Reference - Male"
)
ggplotly(p)Age
tmp2 <- tmp1 %>%
filter(var == "age") %>%
mutate(level = as.factor(level))
p <- ggplot(data = tmp2, aes(x = level, y = OR, colour = country)) +
geom_line(group = 1, position = position_dodge(0.3)) +
geom_errorbar(aes(ymin = `2.5 %`, ymax = `97.5 %`), position = position_dodge(0.3), width = 2) +
geom_point(size = 3, position = position_dodge(0.3)) +
scale_y_continuous(trans = "log10") +
scale_fill_brewer(palette = "Set1")+
labs(title = "Reference - Contunious"
)
ggplotly(p)Education
tmp2 <- tmp1 %>%
filter(var == "education") %>%
mutate(level = as.factor(level))
p <- ggplot(data = tmp2, aes(x = level, y = OR, colour = country)) +
geom_line(group = 1, position = position_dodge(0.3)) +
geom_errorbar(aes(ymin = `2.5 %`, ymax = `97.5 %`), position = position_dodge(0.3), width = 2) +
geom_point(size = 3, position = position_dodge(0.3)) +
scale_y_continuous(trans = "log10") +
scale_fill_brewer(palette = "Set1")+
labs(title = "Reference - No formal education"
)
ggplotly(p)## Warning: Transformation introduced infinite values in continuous y-axis
## Transformation introduced infinite values in continuous y-axis
## Transformation introduced infinite values in continuous y-axis
## Transformation introduced infinite values in continuous y-axis
## Transformation introduced infinite values in continuous y-axis## Warning: `position_dodge()` requires non-overlapping x intervalsMarital
tmp2 <- tmp1 %>%
filter(var == "marital") %>%
mutate(level = as.factor(level))
p <- ggplot(data = tmp2, aes(x = level, y = OR, colour = country)) +
geom_line(group = 1, position = position_dodge(0.3)) +
geom_errorbar(aes(ymin = `2.5 %`, ymax = `97.5 %`), position = position_dodge(0.3), width = 2) +
geom_point(size = 3, position = position_dodge(0.3)) +
scale_y_continuous(trans = "log10") +
scale_fill_brewer(palette = "Set1")+
labs(title = "Reference - Never married"
)
ggplotly(p)## Warning: Transformation introduced infinite values in continuous y-axis
## Transformation introduced infinite values in continuous y-axis
## Transformation introduced infinite values in continuous y-axis
## Transformation introduced infinite values in continuous y-axis
## Transformation introduced infinite values in continuous y-axis## Warning: `position_dodge()` requires non-overlapping x intervalsOccupation
tmp2 <- tmp1 %>%
filter(var == "occupation") %>%
mutate(level = as.factor(level))
p <- ggplot(data = tmp2, aes(x = level, y = OR, colour = country)) +
geom_line(group = 1, position = position_dodge(0.3)) +
geom_errorbar(aes(ymin = `2.5 %`, ymax = `97.5 %`), position = position_dodge(0.3), width = 2) +
geom_point(size = 3, position = position_dodge(0.3)) +
scale_y_continuous(trans = "log10") +
scale_fill_brewer(palette = "Set1") +
labs(title = "Reference - Employed"
)
ggplotly(p)## Warning: Transformation introduced infinite values in continuous y-axis
## Transformation introduced infinite values in continuous y-axis
## Transformation introduced infinite values in continuous y-axis
## Transformation introduced infinite values in continuous y-axis
## Transformation introduced infinite values in continuous y-axis## Warning: `position_dodge()` requires non-overlapping x intervalsLive area
tmp2 <- tmp1 %>%
filter(var == "live_area") %>%
mutate(level = as.factor(level))
p <- ggplot(data = tmp2, aes(x = level, y = OR, colour = country)) +
geom_line(group = 1, position = position_dodge(0.3)) +
geom_errorbar(aes(ymin = `2.5 %`, ymax = `97.5 %`), position = position_dodge(0.3), width = 2) +
geom_point(size = 3, position = position_dodge(0.3)) +
scale_y_continuous(trans = "log10") +
scale_fill_brewer(palette = "Set1")+
labs(title = "Reference - Urban"
)
ggplotly(p)## Warning: Transformation introduced infinite values in continuous y-axis
## Transformation introduced infinite values in continuous y-axis
## Transformation introduced infinite values in continuous y-axis
## Transformation introduced infinite values in continuous y-axis
## Transformation introduced infinite values in continuous y-axis## Warning: `position_dodge()` requires non-overlapping x intervalsMobile
tmp2 <- tmp1 %>%
filter(var == "mobile") %>%
mutate(level = as.factor(level))
p <- ggplot(data = tmp2, aes(x = level, y = OR, colour = country)) +
geom_line(group = 1, position = position_dodge(0.3)) +
geom_errorbar(aes(ymin = `2.5 %`, ymax = `97.5 %`), position = position_dodge(0.3), width = 2) +
geom_point(size = 3, position = position_dodge(0.3)) +
scale_y_continuous(trans = "log10") +
scale_fill_brewer(palette = "Set1") +
labs(title = "Reference - No"
)
ggplotly(p)## Warning: Transformation introduced infinite values in continuous y-axis
## Transformation introduced infinite values in continuous y-axis
## Transformation introduced infinite values in continuous y-axis
## Transformation introduced infinite values in continuous y-axis
## Transformation introduced infinite values in continuous y-axis## Warning: `position_dodge()` requires non-overlapping x intervalsInternet
tmp2 <- tmp1 %>%
filter(var == "internet") %>%
mutate(level = as.factor(level))
p <- ggplot(data = tmp2, aes(x = level, y = OR, colour = country)) +
geom_line(group = 1, position = position_dodge(0.3)) +
geom_errorbar(aes(ymin = `2.5 %`, ymax = `97.5 %`), position = position_dodge(0.3), width = 2) +
geom_point(size = 3, position = position_dodge(0.3)) +
scale_y_continuous(trans = "log10") +
scale_fill_brewer(palette = "Set1") +
labs(title = "Reference - No"
)
ggplotly(p)## Warning: Transformation introduced infinite values in continuous y-axis
## Transformation introduced infinite values in continuous y-axis
## Transformation introduced infinite values in continuous y-axis
## Transformation introduced infinite values in continuous y-axis
## Transformation introduced infinite values in continuous y-axisWe can use anova() function to obtain the P value for each independent variable using likeli-hood ratio test